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Abstract 

Both deterministic and stochastic processes are expected to drive the assemblages of arbuscular mycorrhizal (AIVI) fungi, but 
little is known about the relative importance of these processes during the spreading of toxic plants. Here, the species 
composition and phylogenetic structure of AM fungal communities colonizing the roots of a toxic plant, Ligularia virgaurea, 
and its neighborhood plants, were analyzed in patches with different individual densities of L virgaurea (represents the 
spreading degree). Community compositions of AM fungi in both root systems were changed significantly by the L 
virgaurea spreading, and also these communities fitted the neutral model very well. AM fungal communities in patches with 
absence and presence of L virgaurea were phylogenetically random and clustered, respectively, suggesting that the 
principal ecological process determining AM fungal assemblage shifted from stochastic process to environmental filtering 
when this toxic plant was present. Our results indicate that deterministic and stochastic processes together determine the 
assemblage of AM fungi, but the dominant process would be changed by the spreading of toxic plants, and suggest that 
the spreading of toxic plants in alpine meadow ecosystems might be involving the mycorrhizal symbionts. 
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introduction 

Understanding ecological processes that drive community 
assembly is one of the key ecological topics. Niche-based theory 
emphasizes the importance of deterministic processes (e.g. 
environmental filtering and competitive exclusion) in structuring 
community [1], whereas neutral theory suggests that community 
assembly is mainly determined by the stochastic processes [2]. 
Although the two theories are obviously different, both of them 
have been successfully employed in explaining the diversity and 
assemblage of community. In recent years, more and more 
ecologists have observed that the native communities are 
controlled by both deterministic and stochastic processes, suggest- 
ing that one community composition at any point should be the 
result of interplay of dispersal limitation, environmental filtering 
and species interaction [3,4]. Nonetheless, the current knowledge 
of community assembly is mainly derived from plant communities 
[5,6], whether those unseen communities, such as soil microor- 
ganisms, also follow these assembled rules needs further investi- 
gations. 

Soil microorganisms such as bacteria, fungi and nematodes, are 
very important in regulating biogeochemical cycles [7]. Among 
these microorganisms, arbuscular mycorrhizal (AM) fungi, which 



can form obligately mutualistic associations with the roots of most 
terrestrial plants, are especially renowned because they can 
efficiendy promote plant nutrient uptake [8], improve the stress 
tolerance of their host plants [9] and influence many ecosystem 
processes [10]. Recent evidence has shown that the ecological 
function of an AM fungal community is highly dependent on the 
species composition [11,12], so that exploring the assembly 
mechanisms of AM fungal community is essential before their 
functioning in the fields can be well predicted and even controlled. 
Previous studies have shown that the deterministic factors, such as 
soil nutrient availability [13,14], light inten,sity [15], pH [16], 
phenology [17,18] and host plant identity [19-21], play key roles 
in structuring AM fungal community. However, the stochastic 
processes, such as dispersal limitation [22] and priority effects 
[23,24], can also exert great impacts on the community assembly 
of AM fungi. So far, the assembly mechanisms of AM fungal 
community have been well elucidated using species abundance 
model [23], randomization test [25], neutral model [26] and 
bipartite network analysis [27]. In these studies, however, AM 
fungal species are often considered as phylogenetically equivalent 
entities, and the influence of phylogenetic relatedness between 
species on fungal community assembly is also neglected. In recent 
decades, there has been a growing consensus to merge phyloge- 
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Figure 1. Nonmetric multidimensional scaling (NMDS) ordina- 
tion of AM fungal communities between Ligularia virgaurea 
roots and neighborhood plant roots (a), in L. virgaurea roots (b) 
and in neighborhood plant roots (c) based on Bray-Curtis 
dissimilarity. Plant properties and soil characteristics were fitted as 



vectors onto each ordination plot, and significant vectors at 95% 
confidence level (P<0.05) were displayed onto the NIVIDS ordination 
plots. Control, LD, MD and HD represent no L virgaurea, low density, 
moderate density and high density of L virgaurea, respectively. 
doi:10.1371/journal.pone.0095672.g001 

netic relatedness between species with the studies of community 
ecology, namely the phylogenetic analysis of community structure 
(hereafter referred to as phylogenetic analysis) [28,29]. This 
method can provide great promise for the study of community 
assembly, and demonstrate the importance of phylogenetic 
relatedness between species in driving community assembly. It is 
a key step to clarify whether the functional traits of studied 
organism are conserved (the functional traits are similar among 
closely related species) or not when conducting a phylogenetic 
analysis. Under the assumption of phylogenetic conservatism, in 
theory, the environmental filtering is expected to generate pattern 
of phylogenetical clustering (co-occurring species are more related 
than expected by chance); conversely, competitive exclusion 
should generate phylogenetic overdispersed (co-occurring species 
are more distantly related than expected by chance) [4,30]. The 
stochastic processes should generate a phylogeneticaUy random 
community [4,30]. As the functional traits of AM fungi (e.g. ability 
of hyphal colonization in both roots and soils) are conserved 
[31,32], the phylogenetic analysis might be effective, and by which 
the assembly mechanisms of AM fungal community have been 
partly elucidated by recent studies [31,33,34]. 

Ligularia virgaurea (Asteraceae), a perennial toxic plant for 
livestock, is an indicator species of degraded meadows on the 
Qinghai-Tibetan Plateau [35]. In the past decades, large spreading 
of i. virgaurea in this region has resulted in serious loss of ecological 
and economic services. The large spreading of Z. virgaurea is pardy 
attributed to its aggressive sexual and asexual reproduction [35] 
and exudation of aUelochemicals [36], but whether the spreading 
of L. virgaurea is involving their associated soil microbes remains 
unclear. To fully understand whether the AM fungi have a 
potential role in promoting the L. virgaurea spreading, it is required 
to examine the AM dynamics and explore the ecological processes 
driving the AM fungal assemblage across a spreading gradient, 
because changes of AM fungal assemblage would result in great 
shifts of plant-plant interaction and plant communities [37,38]. 

Recent stadies have suggested that some toxic plants would change 
the interactions between their surrounding plants and belowground 
AM fungi using phytochemicals [39,40], and these changes could 
facilitate the growth of themselves [39,41]. Therefore, we hypothe- 
sized that the spreading of L. virgaurea, like other toxic plants, wotild 
change community assemblage of AM fungi in its neighborhood 
plant roots. To test this hypothesis, we chose four patches with 
different individual densities of L. virgaurea in an alpine meadow 
ecosystem, and the AM fungal communities in the roots of Z. virgaurea 
and its neighborhood plants were analyzed using DNA cloning and 
sequencing techniques. In particular, the objectives of this study were 
to address the foUowdng questions: 1) how do the assemblage of AM 
fungi in the roots of L. virgaurea and its neighborhood plants respond to 
the spreading of L. virgaurerS 2) Which ecological process (environ- 
mental filtering, competitive exclusion or stochastic process) is more 
important in driving the community assembly of AM fungi during the 
spreading of Z. virgaurea? 

Materials and Methods 

Ethics statement 

The field site is supported by Lanzhou University, China. No 
specific permits are required for this study, and also our field works 
did not involve endangered or protected species. 
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Figure 2. The nonlinear curve fitting of AM fungal abundance 
versus species rank in abundance in the roots of Ligularia 
virgaurea (a, ^=9.4, df=6, P=0A5; Kolmogorov-Smirnov 
Z=0.47, P=Q.9S; estimated parameters, 6=1.2, 1 = 7.5) and 
neighborhood plants (b, ^=8.3, df=i^, A'=0.33; Kolmogorov- 
Smirnov Z= 0.40, P= 0.99; estimated parameters, 8 = 9.3, 1 = 3.8). 
Open circles and solid line represent observed values and their fitted 
curve, respectively; filled circles and dotted line represent expected 
values under a neutral model and their fitted curve, respectively. 
doi:1 0.1 371 /journal.pone.0095672.g002 



Filed site and sampling procedure 

This study was conducted in an alpine meadow at the Maqu 
experimental site of the Research Station of Alpine Meadow and 
Wetland Ecosystems of Lanzhou University, which is located in 
the eastern Qinghai-Tibet Plateau of China (33°44'N, 101°52'E; 
3570 m a.s.l.). The climate of this region is cold humid-alpine with 
a mean annual temperature of 1.2°C (11.5°C maximum injuly 
and — 10.0°C minimum in January) and a mean annual 
precipitation of 650 mm. The vegetation in this region consists 
mainly of arctic, alpine and Chinese-Himalayan flora, and is 
dominated with Kobresia capillifolia (Cyperaceae), Elymus nutans 
(Gramineae) and Anemone rivularis (Ranunculaceae). The growing 
season is from May to September. 

Four patches with different individual densities of L. virgaurea, 
including no L. virgaurea (Control), low density (LD), moderate 
density (MD) and high density (HD), were selected within a 
50x50 m flat meadow on 24 August, 2012. These four patches 
could well represent different spreading degrees of L. virgaurea. 
Small and flat research field selected in this study could ensure that 



the original vegetation and soil conditions were relatively 
homogeneous. In each patch, we randomly selected five 
0.5 xO.5 m quadrats. The individual numbers and heights of L. 
virgaurea and the plant species composition in each quadrat were 
measured. In order to determine the AM fungal community 
colonizing L. virgaurea roots in each patch, three individuals of L. 
virgaurea were excavated randomly in each quadrat, and their fine 
roots were collected and pooled as one sample. In each quadrat, 
five sou cores (diameter: 3.8 cm, depth: 25 cm) were taken 
randomly, mixed adequately as one sample in a plastic bag and 
transported to laboratory in an ice box. AH roots excluding L. 
virgaurea were separated carefully from each soil sample and mixed 
thoroughly as one sample to determine AM fungal community 
composition in the roots of neighborhood plants. The L. virgaurea 
roots were distinguished from other plant roots according to their 
morphological character, color and special smell. In order to 
describe easily, all plants in Control patch were also called 
neighborhood plants even there was no L. virgaurea in this patch. 
Each root sample was divided into three subsamples, one for 
determining AM fungal colonization, one for measuring the root 
N and P concentration, and the remaining one for DNA 
extraction. The percent root length colonization by AM fungi 
(%RLC), arbuscular colonization (%AC) and vesicular coloniza- 
tion (%VC) were quantified using the magnified intersection 
method [42] . AH individual plants in each quadrat were clipped to 
the soil surface and grouped into five groups: sedge, grasses, 
legume, forb and L. virgaurea. AU shoots were dried at 80°C for 
48 h, and weighed. All dry shoots of L. virgaurea and non-Z.. 
virgaurea plants in each quadrat were grounded to fine powders for 
determining tissue nutrient concentration, respectively. 

Analysis of soil characteristics and plant N and P 
concentration 

Soil moisture content was measured gravimetricaUy after drying 
at 105°C for 24 h. Soil pH was measured in 1 M KCl (1:5 w/v) 
using a pH electrode. SoU organic carbon and total N were 
analyzed by the CHNS-analysis system (Elementar Analysen 
Systeme, GmbH, Hanau, Germany) with the burning method at 
450°C and 1250°C, respectively. The concentration of soil 
available N (NO3-N+NH4-N) was measured using a FIAstar 
5000 Analyzer (FOSS, Hillerod, Denmark) after extraction with 
2 M KCl (1:5 w/v). Soil available P was extracted following the 
Mehlich-3 method and measured using the molybdate-blue 
colorimetrical method [43]. Plant N content was analyzed by 
the San Continuous Flow Analyzer (Skalar Analytical B.V., 
Breda, Netherlands) using the semi-micro Kjeldahl method. Plant 
P content was measured using the molybdate-blue colorimetric 
method after digestion with sulfuric acid [44] . 

IVlolecular analysis 

DNA was extracted from 35 root samples (20 neighborhood 
plant and 1 5 L. virgaurea root samples, 1 00 mg fine roots for each 
sample) using a Plant DNA Extraction Kit (Tiangen Biotech, 
Beijing, China) following the manufacture's instruction. Partial 
SSU rDNA fragments of AM fungi were amplified via nested PGR 
with a first primer parr GeoA2-Geoll [45] and a second primer 
pair NS31-AML2 [46,47] using the Taq PGR Kit (New England 
Biolabs, MA, USA). The PGR conditions and mixtures of both the 
first and second amplification were the same as the method 
described by Liu et al. [14]. The second PGR products were 
purified using the AxyPrep DNA Gel Extraction Kit (Axygen Inc, 
Hangzhou, China), and the expected DNA products [e. 560 bp) 
were obtained. Thirty-five purified DNA samples were ligated into 
pGEM-T vector (Promega, WI, USA) and transformed into 
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Escherichia coli DH5ot according to manufacture's instructions, 
resulting in 35 clone libraries. For each clone library, 48 putative 
positive transformants were selected randomly and used to prepare 
the plasmid templates using the freezing and thawing method [19]. 
Inserted DNA fragments were re-amplified with primer pair 
NS31-AML2 using the same condition as the second PGR. A total 
of 1572 positive PGR products was obtained from 35 clone 
libraries and then screened using restriction fragment length 
polymorphism (RFLP) with the restriction enzymes HinU and 
//mill (Fermentas, Vilnius, Lithuania). Digested products were 
examined on 2. 5"/!) (w/v) agarose gels with ethidium bromide 
staining, and the RFLP patterns were only compared within one 
sample. One representative clone of each RFLP type in each 
sample was sequenced using vector primer T7 by the Major 
Biotech Gompany (Shanghai, Ghina), for a total of 435 sequences. 
The remaining clones were classified by RFLP patterns and the 
clone numbers of each RFLP/ sequence type were recorded. All 
DNA sequences were edited and compared with public nucleotide 
sequence databases using the online BLAST search tool (http:// 
blast.ncbi.rdm.nih.gov/Blast.cgi), and those with low bit scores and 
high /^-values were suspected as chimeras [48]. Non-AM fungal 
sequences and chimeric sequences were eliminated from the 
dataset before further analysis. All AM fungal sequences were 
multiple ahgned using Glustal W [49], and delimited to species- 
level groups (each species-level group was regarded as an AM 
fungal phylotype) according to 97% sequence similarity using the 
furthest neighbor algorithm in the Mothur program [50]. The 
most closely related sequences from GenBank database and the 
representative sequences of major families of Glomeromycota 
were aligned using the program Glustal W [49], and the neighbor- 
joining tree was constructed using the MEGA 5.0 with the Kimura 
two-parameter model [51]. The classification of our AM fungcd 
phylotypes was according to the phylogeny and taxonomy of 
Glomeromycota (http:/ /schuessler.userweb.mwn.de/amphylo/). 
All representative sequences were submitted to the GenBank 
database under the accession number KF467255-KF467324. 

Statistical analysis 

All statistical analyses were performed using SPSS 13.0 (SPSS 
Inc., Ghicago, IL, USA) or R version 2.15.2 (http://www.r- 
project.org/). AU data were tested for normality after appropriate 
transformation and tested for equality of variance using Levene's 
test. The community composition of AM fungi colonizing roots 
was calculated on the basis of the clone numbers of each phylotype 
in a root sample. To ensure the abundance of dominant and rare 
phylotypes contributed equally to the resultant matrix, the data of 
all AM fungal communities were HeUinger-transformed using the 
function of "decostand" from "vegan" Ubrary of R package [52]. 
Effects of patch or individual density of L. virgaurea on measured 
plant and AM fungal variables were tested by onc-^vay ANOVA, 
and their significant differences were determined using Tukey's 
honestiy significant difference test at the 95"/o confidence level 
(/^0.05). Sampling effort curves of the richness of AM fungal 
phylotypes (Mao Tau) were computed using Estimates 8.0 [53]. 

To elucidate the relationships between AM fungal communities 
and environmental variables, all measured plant and soil variables 
(see Table SI -S3 in File SI) were fitted as vectors onto the non- 
metric multidimensional scaling (NMDS) plots of AM fungal 
communities, which was measured with the Bray-Gurtis dissim- 
ilarity, using the function of 'envfit' from 'vegan' Ubrary of the R 
package [52]. We also tested the correlations between AM fungal 
communities and the matrices of plant species composition, soil 
characteristics, N and P variables using Mantel test, in which the 
matrices of AM fungal communities and environmental datasets 



were represented by the Bray-Gurtis and Euclidean distance, 
respectively [52]. To determine the indicator species of AM fungi 
for each patch, we conducted indicator species analyses (species 
with Indval values ^0.6 are strong indicators) for each AM fungal 
community using the function 'indval' from the 'labdsv' library of 
the R package [54] . 

To test whether the species abundance distribution (SAD) of our 
AM fungal community (community data for L. virgaurea and 
neighboohood plants, respecti\cly) complies with expectation by- 
chance, we fitted Etienne's neutral model using the maximum 
likelihood estimates of metacommunity di\ersity (0) and the 
immigration parameter (I) [55]. Based on the abundance 
distribution of our detected phylotypes, we used the PARI/GP 
algorithm to calculate the G and I, and then run a sequential 
construction scheme that generates 200 expected species abun- 
dance distributions that are consistent with the neutral theory [55]. 
The SAD within the real and simulated data sets was grouped 
using the function "preston" from the R package [56], and the 
goodness-of-fit of our data to the neutral model was assessed using 
test and Kolmogorov-Smimov Z test (in these tests, high P 
values indicate that the data fit the model well). 

To examine the relative importance of spatial effect and 
environmental factors in rc-gulating AM fungal community 
composition, we partitioned the Beta-diversity into spatial and 
environmental components according to the methods described by 
Legendre et al. [57,58]. The spatial structures of the samples were 
modeled using principal coordinates of neighborhood matrix 
(PGNM), which is based on the calculation of a principal 
coordinate analysis (PGoA) of a truncated distance matrix [59]. 
PGNM eigenfunctions represent a spectral decomposition of the 
spatial relationships among samples, and can be considered as 
purely spatial variables in an ordination-based analysis [16]. 
Forward section of PGNM variables was performed, based on 999 
permutation tests (a = 0.05), using the "packfor" library of R 
package [60] . The variation of the community composition data 
was then partitioned between the extracted PGNM spatial 
variables and the soil environmental variables (Table S3 in File 
SI) using the function "varpart" from the "vegan" library of R 
package [52]. The total variations of community compositions 
across four patches were divided into four fractions: (a) variation 
explained by the environmental variables and not spatial 
structured, (b) variation explained by the environmental variables 
and spatial structured, (c) spatially structured variation not 
explained by the environmental variables, and (d) residual 
variation [58]. 

The phylogenetic analysis was performed to explore the 
phylogenetic patterns of our AM fungal communities using the 
"picante" library of R package [61]. The mean nearest taxon 
distance (MNTD) was chosen to quantify the phylogenetic 
relatedness of co-occurring species. We calculated the MNTD 
based on 999 randomizations of the "taxa.labels" null model 
(shuffle distance matrix labels), and the mean nearest taxon 
distances for each species were weighted by species abundance. To 
determine whether the phylogenetic structure of our observed 
community is different from the null community, the nearest taxon 
index (NTI) was calculated based on the differences between the 
obser\ (;d metric and mean metric of the null communities. The 
NTI represents a standardized effect size of MNTD in commu- 
nities, and was calculated using the following formula: 

NTI = - 1 X [MNTDobserved-MNTD,,„dom]/SD(MNTD,^dom), 

where SD represents the standard deviation [31]. The significant 
differences between NTI and null expectation of zero were tested 
using two-tailed T test at 95% confidence level: if the NTI is 
significantiy higher than zero, this indicates that the community 



PLOS ONE I www.plosone.org 



4 



April 2014 | Volume 9 | Issue 4 | e95672 



Ecological Processes Driving AM Fungal Assemblage 



Table 1. Relative abundance (%, proportion of clone numbers) of each arbuscular mycorrhizal fungal phylotype and indicator 
species in the roots of Ligularia virgaurea and its neighborhood plants in different patches. 
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Control, LD, MD and HD represent no 1. virgaurea, low density, moderate density and high density of L virgaurea,, respectively; u: uncultured. 
^The GenBank accession number of the representative sequence. 

"^'new' in parentheses indicates <97% similarity with published sequences in the GenBank database. 

*^Bold font indicates that the phylotype Is a significant (lndval>:0.6, P<0.05) indicator species for a particular patch. 

doi:1 0.1 371 /journal.pone.0095672.t001 



phylogenetic structure should be clustered; if the NTI is lower than 
zero, the community should be phylogenetically overdispersed; 
and if the NTI does not differ from zero, the community should be 
phylogenetically random. 

Results 

Plant and soil properties in different patches 

The individual densities (one-way ANOVA, F= 174.7, 
P<0.001) and heights (i^= 5.2, P=0.02) of L. virgaurea were 
gradually increased from LD to HD patches (Table SI in File SI), 
suggesting that our selected patches could well represent the 
spreading degree of L. virgaurea. With the spreading of L. virgaurea, 
both species richness [F— 6.3, P= 0.005) and shoot biomass 
(^^=5.1, P= 0.012) of their neighborhood plants were declined 



significandy (Table SI in File SI), and these changes were 
especially remarkable for the forb and grass (Table S2 in File SI). 
Plant tissue P concentration in both L. virgaurea and their neighbors 
varied among patches, with the highest P concentration synchro- 
nously occurred in HD patch (Table SI in File SI). The lowest soil 
available N and N/P ratio were both present in HD patch, 
whereas other soil characters (soil moisture, pH, soil C/N ratio, 
etc.) were similar across the density gradient of L. virgaurea (Table 
S3 in File SI). 

AM fungal colonization and community composition 

The spreading of L. virgaurea significandy changed %AC 
(/'=3.9, P=0.03) and %VC (i='=6.3, P=0.01) in their neighbor- 
hood plant roots, but did not affect the AM colonization in the 
roots of themselves (Table S4 in File SI). A total of 324 AM fungal 
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Table 2. Mantel tests between matrices of AM fungal community in the roots of Ligularia virgaurea and its neighborhood plants 
and matrices of N and P variables, plant community and soil characteristics. 





Matrices 


Component variables 


Ligularia virgaurea roots 


Neighborhood plant roots 






Mantel r 


lvalue 


Mantel r 


lvalue 


Plant community 


Species compositions 


0.21 


0.049 


0.24 


0.009 


P availability 


Soil AP, shoot P, root P 


0.18 


0.11 


0.23 


0.019 


N availability 


Soil AN, shoot N, root N 


0.08 


0.263 


0.03 


0.371 


Soil characteristics 


All measured soil variables 


0.18 


0.084 


-0.04 


0.682 



doi:10.1371/journal.pone.0095672.t002 



sequences (74.5% of total sequences), based on 1184 AM fungal 
clones, were obtained in this study. These sequences could be 
delimited into 30 phylotypes with sequence identity £97% 
(Table 1, Figure SI in File SI). Of these phylotypes, seven were 
related to sequences of described species, six to spore-derived 
sequences, fifteen to uncultured AM fungi, and the remaining two 
phylotypes (Glo-14 and Aca-17) were undescribed previously 
(<97% identity with pubhshed sequences; Figure SI in File SI). 
These phylotypes belonged to 10 Glomeromycotan families: 
fourteen Glomeraceae, four Acaulosporaceae, four Claroideoglo- 
meraceae, two Archaeosporaceae, one Gigasporaceae, one 
Ambisporaceae, one Pacisporaceae, one Diversisporaceae, one 
Geosiphonaceae and one Paraglomeraceae phylotype (Figure SI 
in File SI). Sampling effort curves clearly showed that the majority 
of AM fungal phylotypes in both root systems were sufficiently 
identified (Figxire S2 in File SI). In total, 26 and 23 AM fungal 
phylotypes were detected in L. virgaurea and neighborhood plant 
roots, respectively, and 18 phylotypes were same in both root 
systems (Table 1). In addition, the community compositions of AM 
fungi in both root systems were similar (Figure la) and correlated 
with each other significantly (Mantel r=0.57, P= 0.001). 

Spreading ot L. virgaurea did not change the AM fungal richness 
and Shannon- Weiner diversity in the roots of L. virgaurea (richness: 
F= 0.5, P= 0.62; diversity: F= 0.9, P= 0.43) and neighborhood 
plants (richness: P=0.79; diversity: i?=0.3, P=0.80). 

However, the community compositions of AM fungi in both root 
systems varied among patches (Figure lb,c). An average of 42% 
(neighborhood plant roots) and 48% (L. virgaurea roots) abundance 
per sample belonged to one phylotype, and the most dominant 
phylotypes in each root system-patch combination were different, 
with exception of the phylotype Glo-2, which was predominant in 
neighborhood plant roots in both MD and HD patches (Table 1). 
In addition, some AM fungal phylotypes were only detected in a 
particular patch or root system-patch combination; for example, 
the Glo-1 was specific for HD patch, and Scu-21 was only detected 
in L. virgaurea roots in HD patch (Table 1). Indicator species were 
also different between AM fungal communities in two root 
systems; for example, Glo-10 was the indicator species for 
neighborhood plants in HD, but Glo-9 for L. virgaurea in MD 
(Table 1). 

The community composition of AM fungi in L. virgaurea roots 
was significantly related with L. virgaurea density (r = 0.6, 
P= 0.004), soil available N/P ratio (/ = 0.57, P=0.01) and soil 
C/N ratio (r^ = 0.43, P—Q.Q'i); of which, the L. virgaurea density 
was the most relevant variable (Figure lb). Some soil and plant 
variables were also correlated with the AM fungal communities 
colonizing neighborhood plant roots (Figure Ic), and the root N/P 
ratio (/ = 0.45, P= 0.003) was the most relevant variable, followed 
by root P (/=0.39, P=0.01) and L. virgaurea density (/ = 0.37, 



P= 0.02). Mantel test revealed that the plant species composition 
was significandy correlated with AM fungal communities in both 
root systems (Table 2). P availability was only related with AM 
fungal community in neighborhood plant roots, while N 
availability and soil characteristics were not correlated with both 
fungal communities (Table 2). 

Neutral model fitting, variation partitioning and 
phylogenetic analysis of AM fungal community 

We plotted the association between AM fungal relative 
abundance and species rank in abundance as predicted by a 
neutral model to test stochastic process in driving community 
assembly (Figure 2). The abundance distribution of our AM fungal 
phylotypes in the roots of Z. virgaurea 9 A, df=6, P=0.15; 
Kolmogorov-Smirnov ^= 0.47, P= 0.98; estimated parameters, 
9 = 1 .2, 1 = 7.5) and neighborhood plants (/ = 8.3, df=4,P= 0.33; 
Kolmogorov-Smirnov ^= 0.40, P= 0.99; estimated parameters, 
9 = 9.3, I = 3.8) showed a close fitting to the neutral model. 

Analysis of variation partitioning explained 21.8% of the 
variation of AM fungal community composition, based on 
adjusted value (adj. P^), in the roots of L. virgaurea. Of these 
variation, 17.5% (adj. P^a+b = 0.175) was attributed to environ- 
mental variables, and 4.4% (adj. P^f^-c = 0.044) to spatially 
structured variation; moreover, environmental and spatial vari- 
ables jointiy explained 17.3% of the variation of AM fungal 
community composition in the roots of neighborhood plants; of 
which, 16.9% (adj. P^,+,, = 0.169) and 9.1% (adj. P\+, = 0.091) 
were attributed to environmental and spatial factors, respectively. 

Our phylogenetic analysis successfully revealed the phylogenetic 
patterns of AM fungal community in each patch. The NTI in 
three patches (LD, MD and HD) was significantly higher than zero 
in both neighborhood plant roots and L. virgaurea roots, whereas it 
was not significantly different from zero in Control patch 
(Figure 3). These results suggest that AM fungal communities in 
LD, MD and HD patches were phylogenetically clustered, and the 
community in Control patch was phylogenetically random. 
Although not statistically significant, the NTI for both root 
systems declined gradually from LD to HD (Figure 3). 

Discussion 

In the present study, we did not observe significant changes of 
AM fungal richness and diversity during the spreading of L. 
virgaurea. However, a total of seven AM fungal phylotypes were 
only present in HD patch, suggesting that these fungi might be 
immigrants. It is possible that these immigrated fungi might move 
from deep soil layer [62] or be introduced from the excrements of 
soil arthropods [63]. This speculation is also supported by our 
analysis of AM fungal communities: the community composition 
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Ligularia virgaurea Neighborhood Plants 



Figure 3. Nearest taxon index (NTI, mean ± SE) of arbuscular mycorrhizal fungal community in the roots of Ligularia virgaurea and 
its neighborhood plants in different patches. Mean NTI is significantly higher than zero, indicating that the community is more 
phylogenetically clustered than expected by chance. Asterisk indicates that the NTI has significant differences with zero at 95% confidence level 
(PsO.05). Control, LD, MD and HD represent no L virgaurea, low density, moderate density and high density of L virgaurea, respectively. 
doi:1 0.1 371 /journal.pone.0095672.g003 



in HD was especially distinct with that in Control, LD and MD 
patches, whilst there was less difference among the later three 
patches. These findings partly support our hypothesis, and are in 
line with previous studies, in which some forbs could change the 
community composition of AM fongi colonizing their surrounding 
plants [21,64]. In addition, significant changes of most soil and 
plant variables only occurred in the HD but not in other patches. 
This might be another explanation for the distant community 
composition between HD with others, and indicates that the AJvI 
fungal community would be largely changed when the individual 
density of L. virgaurea exceeds a threshold value. 

Our Mantel test revealed that the community composition of 
AM fungi was highly correlated with aboveground plant 
community, highlighting the importance of plant species compo- 
sition in regulating their associated fungal community. Further- 
more, some soil characteristics (e.g. soil C/N ratio and P 
availability) were also related with AM fungal communities in 
both L. virgaurea and neighborhood plant roots, suggesting that 
abiotic factors could also exert effects on the mycorrhizal 
symbionts. These findings agree with many previous studies 
showing that the AM fungal communities were determined by 
both plant communities [65] and some environmental factors 
[13,14], and highlight the importance of deterministic process 
based on niche differentiation in driving the community assembly 
of AM fungi. In fact, changes of plant community and soil 
properties in this study should be ultimately attributed to the 
spreading of Z. virgaurea, because its spreading could suppress the 
seedling recruitment of neighborhood plants by allelochemicals 
[36] as well as change soil characteristics via changing the activities 
of soil enzymes [66] . Further studies are encouraged to address 
how L. virgaurea spreading directiy and indirecdy regulates the AM 
fungal assemblages in alpine meadow ecosystems. 

Our AM fungal communities were typically over-dominated by 
a single phylotype, and these dominant phylotypes varied in most 
patch-root system combinations. Similar results were also found in 
other studies, showing the importance of priority effects in fungal 



"meta-community" assembly [23,24]. It has been suggested that 
the overdominance of AM fungi in a community might result from 
a positive feedback occurring during the build-up of the plant- AM 
fungal community [23]. That is, a 'founder AM fungal species' 
colonized plant roots during earlier community succession would 
benefit from more plant derived-carbon than "latecomers", which 
would favor its growth and spreading through the soil, and 
increase its probability of colonizing newly formed roots. Our 
neutral model fitting also confirms the role of stochastic process in 
structuring AM fungal community. Collectively, our results 
corroborate other findings [16,67,68], and suggest that determin- 
istic and stochastic processes jointly determine the community 
assembly of AM fungi. In our case, however, a stronger effect of 
environmental component on the community composition was 
demonstrated by the variation partitioning, indicating that the 
effect of deterministic process was stronger than stochastic process 
in structuring AM fungal community during the spreading of L. 
virgaurea. 

It is interesting that the AM fungal communities in patches with 
absence (Control) and presence (LD to HD) of L. virgaurea were 
phylogenetically random and clustered, respectively. These results 
suggest that the community assembly of AM fungi in patch 
without L. virgaurea was mainly determined by stochastic process, 
but once this toxic plant was present, the environmental filtering 
became as a major process in structuring fungal communities. 
Competitive exclusion can also produce phylogenetically clustered 
pattern even when the functional traits is evolutionally conserved 
[69] , but in our case, the shift of phylogenetic pattern should be 
attributed to environmental filtering from the allelochemicals 
secreted by L. virgaurea, because the major difference between 
Control and other patches was the presence or absence of L. 
virgaurea. Although we did not measure the effects of i. virguarea's 
allelochemicals on fungal community, previous studies have shown 
that the allelochemicals could significandy influence fungal 
community composition by suppressing the spore germination of 
some AM fungal species [70] and the formation of mycorrhizal 
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association [39,40] . It is therefore that the allelochemicals of L. 
virgaurea in soil might be an environmental stress for mycorrhizal 
fungi, and this can in part explain the phylogenetically random 
pattern observed in our Control patch, where the AM fungi was 
free from allelochemicals of L. virgaurea. Nonetheless, further 
researches are required to address how the allelochemicals of L. 
virgaurea affect the assemblage and functioning of AM fungal 
community, and only by which can we fully understand whether 
the spreading of L. virgaurea is involving the indirect effects of 
cdlelochemicals on mycorrhizal communities of their surrounding 
plants. 

In conclusion, our study shows that the spreading of L. virgaurea 
in an alpine meadow ecosystem could significantly change the 
assemblage of AM fungi inside their surrounding plant roots. The 
phylogenetic structure of AM fungal community dramatically 
changed from phylogenetical random to clustering during the 
spreading of L. virgaurea. The community assembly of AM fungi in 
patch without L. virgaurea was mainly determined by stochastic 
process, whereas the environmental filtering became as a major 
process once the L. virgaurea was present. Our findings indicate that 
the spreading of L. virgaurea in alpine meadow ecosystems might 
involve with the mycorrhizal symbionts, and also highlight the 
importance of phylogenetic relatedness between species in driving 
AM fungal community composition. Further studies are encour- 
aged to address the functional changes of AM fungal communities 
during L. virgaurea spreading, as well as the feedbacks of changed 
AM fungal communities on aboveground plant communities. 

Supporting Information 

File SI Contains Table SI, Changes in plant properties for 
Ugularia virgaurea and its neighborhood plant community in 
different patches. Data are means ± SE (n = 5). Significant 
differences among patches within each variable were tested using 
Tukey's honestiy significant difiFerence test (-PsO.05) and are 
indicated by dissimilar letters. Table 82, Changes in plant richness 
and shoot biomass in different patches. Data are means ± SE 
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